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Viscosity and magnetic fields drive diflerentially rotating stars toward uniform rotation, and this 
process has important consequences in many astrophysical contexts. For example, merging binary 
neutron stars can form a "hypermassive" remnant, i.e. a diff'erentially rotating star with a mass 
greater than would be possible for a uniformly rotating star. The removal of the centrifugal support 
provided by differential rotation can lead to delayed collapse of the remnant to a black hole, accom- 
panied by a delayed burst of gravitational radiation. Both magnetic fields and viscosity alter the 
structure of differentially rotating stars on secular timescales, and tracking this evolution presents 
a strenuous challenge to numerical hydrodynamic codes. Here, we present the first evolutions of 
rapidly rotating stars with shear viscosity in full general relativity. We self-consistently include vis- 
cosity in our relativistic hydrodynamic code by solving the fully relativistic Navier-Stokes equations. 
We perform these calculations both in axisymmetry and in full 3-1-1 dimensions. In axisymmetry, the 
resulting reduction in computational costs allows us to follow secular evolution with high resolution 
over dozens of rotation periods (thousands of M). We find that viscosity operating in a hyper- 
massive star generically leads to the formation of a compact, uniformly rotating core surrounded 
by a low-density disk. These uniformly rotating cores are often unstable to gravitational collapse. 
We follow the collapse in such cases and determine the mass and the spin of the final black hole 
and ambient disk. However, viscous braking of difi^erential rotation in hypermassive neutron stars 
does not always lead to catastrophic collapse, especially when viscous heating is substantial. The 
stabilizing infiuences of viscous heating, which generates enhanced thermal pressure, and centrifu- 
gal support prevent collapse in some cases, at least until the star cools. In all cases studied, the 
rest mass of the resulting disk is found to be 10-20% of the original star, whether surrounding a 
uniformly rotating core or a rotating black hole. This study represents an important step toward 
understanding secular effects in relativistic stars and foreshadows more detailed, future simulations, 
including those involving magnetic fields. 

PACS numbers: 04.25.Dm, 04.40.Dg, 97.60.Jd 



I. INTRODUCTION 

The field of numerical relativity has matured to a stage 
where it is possible to simulate realistic systems of astro- 
physical interest. In this paper, we examine the global 
effects of viscosity on differentially rotating, relativistic 
stars. Viscosity can have significant effects on the sta- 
bility of neutron stars. For example, it can drive a sec- 
ular bar instability in rapidly rotating neutron stars, as 
shown in Newtonian gravitation P, and in general rel- 
ativity jSj . Viscosity can suppress the r- modes 0, 1^ and 
other gravitational-radiation driven instabilities, includ- 
ing the secular bar modes . Viscosity also destroys dif- 
ferential rotation, and this can cause significant changes 
in the structure and evolution of differentially rotating 
massive neutron stars. 

Differentially rotating neutron stars can support sig- 
nificantly more rest mass than their nonrotating or uni- 
formly rotating counterparts, making "hypermassive" 
neutron stars possible 0, @|- Such hypermassive neu- 
tron stars can form from the coalescence of neutron star 
binaries 0, ITol [lTI| or from rotating core collapse. The 
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stabilization arising from differential rotation, although 
expected to last for many dynamical timescales, will ul- 
timately be destroyed by magnetic braking and/or vis- 
cosity These processes drive the star to uniform 
rotation, which cannot support the full mass of the hy- 
permassive remnant. This process can lead to "delayed" 
catastrophic collapse to a black hole, possibly accompa- 
nied by some mass loss. Such a delayed collapse might 
emit a delayed gravitational wave signal detectable by 
laser interferometers. Moreover, the collapse, together 
with any residual gas in an ambient accretion disk, could 
be the origin of a gamma-ray burst (GRB). 

Both magnetic fields and viscosity can destroy differen- 
tial rotation in a rapidly rotating star 12, 13, 14]. Sim- 
ple estimates show that the magnetic braking (Alfven) 
timescale for a laminar field is much shorter than the 
timescale of molecular (neutron) viscosity in a typical 
massive neutron star. Hence magnetic fields are expected 
to be the principal mechanism driving neutron stars to- 
ward rigid rotation. Phase mixing arising from magnetic 
braking [l^llSll. or other possible magnetohydrodynamic 
instabilities |l5l Il6l | might stir up turbulence. Turbu- 
lent shear viscosity could then dominate the subsequent 
evolution. In this paper, we arc primarily interested in 
identifying the global evolutionary consequences of shear 
viscosity in a relativistic star, independent of the detailed 
nature or origin of the viscosity. 
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To explore the consequences of the loss of differential 
rotation in equilibrium stars, we study the secular evolu- 
tion of differentially rotating relativistic stars in the pres- 
ence of a shear viscosity. Viscosity and magnetic fields 
have two things in common: (1) they both change the 
angular velocity profiles of a differentially rotating star, 
and (2) they both act on secular timcscales, which can 
be many rotation periods. The latter inequality poses a 
severe challenge to numerical simulations using a hydro- 
dynamic code. It is too taxing for a hydrodynamic code 
using an explicit differencing scheme to evolve a star for 
physical realistic secular timescales. To solve this prob- 
lem, we artificially amplify the strength of viscosity so 
that the viscous timescale is short enough for numeri- 
cal treatment. However, we keep the viscous timescale 
substantially longer than the dynamical timescale of the 
stars, so that the evolution of the star remains quasi- 
stationary. We then check the validity of our results 
by reducing the viscosity on successive runs and test- 
ing that the viscosity-induced physical behavior is un- 
changed; rather, only the timescale changes and does so 
inversely with the strength of viscosity. A more detailed 
discussion of the expected scaling is presented in Sec- 
tion |IIE1 (IJ] . 

To study viscous evolution, we need to perform long 
simulations in full general relativity. Typically, we evolve 
the stars in axisymmetry. This allows us to follow the 
secular evolution of the stars with high resolution in a 
reasonable amount of time. Viscosity can, however, drive 
nonaxisymmetric instabilities when a star is rapidly ro- 
tating. To test for such instabilities, we also perform 
lower-resolution, three-dimensional (3D) simulations on 
the most rapidly rotating stars we consider. 

For non-hypermassive neutron stars that are slowly 
and differentially rotating, we find that viscosity sim- 
ply drives the whole star to rigid rotation. If the non- 
hypermassive neutron star is rapidly and differentially 
rotating, however, viscosity drives the inner core to rigid 
rotation and, at the same time, expels the material in 
the outer layers. The final system in this case consists 
of a rigidly-rotating core surrounded by a low-density, 
ambient disk in quasi-stationary equilibrium. 

Our most interesting results concern the fate of hy- 
permassive neutron stars. We numerically evolve four 
models with different masses and angular momenta. We 
find that in all cases, viscosity drives the cores to rigid ro- 
tation and transports angular momentum outwards into 
the envelope. As a result, the core contracts in a quasi- 
stationary manner, and the outer layers expand to form 
a differentially rotating torus. Of the four models we 
have studied, the star with the highest mass collapses to 
a black hole, with about 20% of the rest mass leftover to 
form a massive accretion disk. On the contrary, the other 
three stars do not collapse to black holes, but form star -|- 
disk systems, similar to the final state of the rapidly ro- 
tating non-hypermassive neutron stars described above. 
As will be discussed in Section III Fl viscosity generates 
heat so that the stars do not evolve adiabatically in gen- 



eral. The extra thermal pressure due to viscous heating 
helps to support the stars. We also consider the hmit of 
rapid cooling, whereby the heat generated by viscosity is 
immediately removed from the stars. Of the three stars 
which do not collapse to black holes in the no-cooling 
limit, we found that the one with the lowest angular 
momentum undergoes catastrophic collapse in the rapid- 
cooling limit. About 10% of the rest mass is leftover to 
form an accretion disk in this case. To test the validity 
of the axisymmetric results, we perform 3D simulations 
to check for any nonaxisymmetric instabilities. We do 
not find any unstable nonaxisymmetric modes and the 
3D results agree with the axisymmetric results. 

Our results suggest that viscous braking of differential 
rotation in a hypermassive neutron star can, but does 
not always, lead to catastrophic collapse. When catas- 
trophic collapse does occur, the remnant is a black hole 
surrounded by a massive accretion disk. This outcome 
is very different from that of the collapse of an unstable, 
rigidly-rotating "supramassive" neutron star, in which 
the whole star collapses to a black hole, leaving only a 
tiny amount of material to form a disk 0, ^| . Many 
models for GRBs require a massive disk around a rotat- 
ing black hole to supply energy by neutrino processes |22j . 
Our results suggest that viscous forces in a hypermassive 
star could lead to the formation of a massive disk around 
such a black hole. 

The structure of this paper is as follows. In Section ITU 
we derive the relativistic Navier-Stokes equations con- 
taining shear viscosity in a 3-1-1 form suitable for numer- 
ical integration, and describe how we evolve them in both 
axisymmetry and full 3+1 dimensions. We then describe 
in Section III several tests that we perform to check our 
code. We present the results of our simulations on five se- 
lected stars in Section IV. Finally, we briefly summarize 
and discuss our conclusions in Section V. 



II. FORMALISM AND NUMERICAL METHODS 



A. Evolution of the gravitational fields 

Throughout this paper, Latin indices denote spatial 
components (1-3) and Greek indices denote spacetime 
components (0-3). We adopt geometrized units, so that 
G = c = 1. We evolve the 3- metric 7ij and the extrin- 
sic curvature Kij using the BSSN formulation The 
fundamental variables for BSSN evolution are 



ln[det(7i 



lij = 

K EE 

Aij = 



12 



(1) 

(2) 
(3) 

(4) 

(5) 
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The evolution and constraint equations for these fields 
are summarized in |24| (hereafter Paper I). In the pres- 
ence of matter, these evolution equations contain the fol- 
lowing source terms: 



P = 

Si-i — 



(6) 



where T"'^ is the stress tensor, and ria — (—a, 0,0,0) is 
the future-directed unit normal to the time slice. One 
must impose gauge conditions which specify the lapse a 
and the shift We use a K-driver lapse and Gamma- 
driver shift, as described in Paper I. The numerical imple- 
mentation of the equations is discussed in Paper I, with 
some improvements to enhance stability described in 25]. 
The latter are particularly relevant for the post-collapse 
version of our code that we implement with black-hole 



give the relativistic continuity, energy, and Navier-Stokes 
equations 



dtp* 
dte* 
dtSk 



+ d,{p*v') - 



(15) 



+ d,{e*v') = -ae«^(poe)('-^^/^a"'5a„^ (16) 



1 



(17) 



where = jvP is the 3- velocity. The quantity vP is 
determined by the normalization condition u'^u^, — — 1, 
which yields 



= pl + e-^^f^SS, 



r-i 



, (18) 



B. 3+1 relativistic Navier-Stokes equations 

We treat the matter in our neutron stars as an imper- 
fect fluid with a shear viscosity, but no bulk viscosity and 
no heat conduction. The stress tensor for the fluid is 



{pa + Poe + P)u^u^ + Pg^^ - 2ria^„ . (7) 



Here, poj and are the rest-mass density, specific 
internal energy, pressure, and fluid four-velocity, respec- 
tively. The quantity rj is the coefficient of viscosity and 
is related to the kinematic viscosity by 77 = pqv. The 
shear tensor (7^,^ is defined by 

C^tiv = + a{tiUy) - -U°';aigfj.v + U^U„) , (8) 

where is the fluid 4-acceleration. We assume a F-law 
equation of state 



P = (F - l)poe . 
Our fundamental fluid variables are 



e* = {poef'^au^e^'>' 
Sk = p*huk , 



(9) 



(10) 

(11) 
(12) 



where h — 1 + e + P/po is the specific enthalpy. The 
conservation of stress-energy 

T^"'.^ = (13) 

and the law of baryon number conservation 

V^(pow^) = (14) 



where w = pi^avP . 

The stress-tensor T^'^ generates the following source 
terms in the field evolution equations: 



= hwe~^'>'-P 
2n 



''Si - — {(Tti - CTijf]') , 

a 



(19) 

(20) 
(21) 



C. 2+1 relativistic Navier-Stokes equations 



Many of the systems we evolve possess and maintain 
symmetry about their rotation axis, which we set to be 
the z-axis. Then we can eliminate one dimension and 
simplify the equations. We utilize axisymmetry and fol- 
low |27l l2q | to evolve the field and hydrodynamic vari- 
ables on the 2/ = plane. The data off of this plane can 
be obtained by rotating the data on this plane. As we 
explain in Section III G II we find it advantageous when 
performing 2+1 simulations to evolve the hydrodynamic 
equations H15|l - H17|) in cylindrical coordinates but on a 
Cartesian {xz) grid. On the y = Q plane, the cylindrical 
coordinates m = + J/^, ^, and f — arctan(y/a;) are 
related to the Cartesian coordinates cc, y, and z as fol- 
lows: Tu ^ X, wLp ^ y, z ^ z. Using these relations. 
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Eqs. (|15|I - H17|) in cylindrical coordinates can be written 



(22) 



dt{SA - 2ae^'^T](j\) + -dsixSAV^ - 2xae^''' r^a^ 



1 ,^ 



-{SyV^ - 2xae'"^'qal)5Ax - ae'"^dAP (23) 



+ae''^g"f^A 



-pohuaUp + riaai3 



dt{Sy-2ae^U^l) 



(24) 



(25) 



^^dsix'SyV^ - 2ae^*r/:EV^) = 

X 



where the indices A and B run over x and z (c.f. 
Eq.(2.10)-(2.13) of ;28]). 



D. Hierarchy of timescales 

There are two dynamical timescales for a rotating star. 
Gravity provides the free-fall timescale rpp 



rpF 



M 



\ 0.2 



) \2Mc. 



s , (26) 



where M is the gravitational mass of the star (or merged 
binary remnant) and R is the radius. If the star is rotat- 
ing, its rotation period Prot provides another important 
timescale: 



P.ot = ^ = 1.6 X 10- 



4000s- 



(27) 



Dynamical instabilities (e.g. instability to radial collapse 
or to dynamical bar formation) will act on the above 
timescales. 

The stars we study are dynamically stable initially, so 
their structure is altered on secular timescales. Rotating 
compact stars may be secularly unstable to gravitational- 
radiation driven instabilities. The strongest instabilities 
of this kind are the (nonaxisymmetric) r-modes and the 
bar mode. The timescale of the I — m — 2 r-mode 
instability is given by 



50 



UoOOs-i J V 0-2 j V2M0 



s . (28) 



The gravitational-radiation driven (Dedekind) bar-mode 
instability occurs if the star is rapidly rotating so that 
T/|Ty| > /3s, where is the ratio of kinetic to grav- 

itational potential energy . The threshold (3g sa 0.14 
for a Newtonian star (M/R <C 1), and decreases as the 



compactness of the star (i.e. M/R) increases *2^. The 
timescale of this instability is estimated to be _30J 



GW 
bar 



V 0.2 ) \2Mq 



fT/\W\ 



f3s 



0.1 



(29) 

Viscosity alone can also drive a (Jacobi) bar-mode in- 
stability. The threshold is identical for a Newtonian star 
{Ps ~ 0.14) but increases as the compaction increases 
The relevant timescale is (see [H, p. 99) 



bar 



R^ f T/\W\-f3s 



0.1 



(30) 



R^/v 



This is comparable to the viscous timescale 
discussed below. 

Magnetic fields coupled to the matter will redistribute 
angular momentum. In fact, even an initially small mag- 
netic field frozen into the matter will be wound up and 
can destroy differential rotation in the star on an Alfven 
timescale 



yy ameren 



100 



B 



1012G 



' (m/rV'^ 



(31) 



Viscosity will also redistribute angular momentum on 
a viscous timescale Tvis- One form of viscosity present 
in neutron stars is due to the transport of energy 
and momentum of neutrons. This viscosity acts on a 
timescale I3lll33 



10**Tq 



/ M 



V 2Mp 



9/2 



-23/4 



(32) 



where Tg = T/IO'^K , and T is the characteristic temper- 
ature. 

It is widely believed that, for cold neutron stars (T < 
lO^-ftT), the neutron fluid in the inner crust condenses into 
a superfluid of Cooper pairs j33j |. while in the interior, 
the neutrons could form a.^P2 superfluid (although 
this is less certain), and the protons a "'^So superfluid. In 
the case of neutron superfluidity, r„_vis will vanish, and 
the dominant viscosity will be due to electron-electron 
scattering js^H^ 



IO^Tq^ 



M 
2Mp 



s . 



(33) 



Electron and proton fluids are forced to move together 
in the MHD limit ^36,]. Differences in velocity between 
the neutron and proton-electron fluids are damped fairly 
quickly by mutual friction [s^ Is^ . 

Viscosity can be used as a model for turbulence in cer- 
tain situations. Turbulence may occur in young neutron 
stars as a result of pure hydrodynamic effects or magnetic 
instabilities . Turbulence is often modeled by the "a- 



disk" law, in which a shear stress T, 



-aP is added 

to the hydrodynamic equation (see e.g. [3^, Chap. 14). 
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Here a is a nondimensional constant (which should not 
be confused with the lapse function) with values in the 
range 0.001 < a < 1. The viscosity in this model roughly 
corresponds to ~ ^turbi'turb ~ aRcs, where Wturb is the 
velocity of turbulent cells relative to the mean fluid mo- 
tion, Zturb is the size of the largest turbulent cell, and Cs 
is the sound speed. The corresponding timescalc is 

- -PF - ^ f^y''7^) s , (34) 
a a V 0.2 y \2Mq) ' ^ ' 

which is much shorter than all the other secular 
timescales. Hence turbulent viscosity, if present, is likely 
to dominate the secular evolution of differentially rotat- 
ing stars. 

Thermal energy is radiated away primarily by neutri- 
nos. For hot neutron stars (T > lO^-fiT), the cooling is 
dominated by the direct URCA process, and the star 
cools on a timescale TcooI ^ lO^Tg"'* s (see js^ and [ssj . 
Chap. 11). For cooler neutron stars, the cooling is domi- 
nated by the modified URCA process, and the star cools 
on a timescale TcooI ~ lO^Tg"^ s Depending on the 
temperature and the nature of the viscosity, the cooling 
timescale may be greater than or less than the viscous 
timescale. If TvIs <C TcooI, then the heat generated by 
viscosity will build up inside the star. Otherwise, it will 
be radiated away as quickly as it is generated. We study 
both limits in this paper. 



E. Dynamically modeling secular effects 

Secular effects will in general take many rotation peri- 
ods to significantly affect the structure or velocity profile 
of a differentially rotating star. This poses a challenge to 
the numerical treatment of these changes. Because of the 
short Courant timestep required for numerical stability, 
it would be computationally prohibitive to evolve a star 
for such a long time using an explicit finite differencing 
scheme. The use of a fully implicit scheme for the finite 
differencing can allow stable evolutions with larger AT. 
Each timestep is, however, much more computationally 
expensive as it involves matrix inversion. Moreover, no 
fully implicit routine for the coupled Einstein field and 
relativistic hydrodynamic equations exists at present. 

The secular timescales are so much longer than the 
dynamical timescales that the star can be thought of as 
evolving quasi-statically. Therefore, it might be possible 
to treat the secular evolution in the quasistatic approxi- 
mation, as in typical stellar evolution (Henyey) codes, by 
constructing a sequence of equilibrium configurations up 
to the moment that stable equilibrium can no longer be 
sustained. This approach has been used to study the vis- 
cous evolution of differentially rotating white dwarfs ^3 ■ 
However, building the required equilibrium models in full 
general relativity is a nontrivial task. It would be par- 
ticularly difficult to identify the meridional currents and 
core-halo bifurcation that often arise in rapidly rotating 



configurations. More significantly, it would not be pos- 
sible to follow the evolution of the configuration with a 
quasistationary approach if a dynamical instability (i.e. 
collapse) is triggered during the secular evolution. 

Instead, we use our relativistic hydrodynamic code 
and artificially amplify the strength of viscosity so that 
the viscous timescale is short enough to make numer- 
ical treatment tractable. However, we keep the viscous 
timescale sufficiently long that the hierarchy of timescales 
is maintained, and the secular evolution still proceeds in 
a quasi-stationary manner. The behavior of the real sys- 
tem can then be determined by rescaling the time vari- 
able to adjust the viscous timescale to its physical value. 
The characteristic viscous timescale is 

Tvis-pi?' <r;>-i , (35) 

where < > is an averaged value of r/ across the star. 
Suppose we evolve the same star, once with t^is — 
and once with Tvis — T2. If both ri and T2 are large 
enough so that they do not compete with the dynamical 
timescale, but shorter than any other secular timescale 
(see, e.g. Section lllD() . then the configuration of the star 
with viscosity ri at time t will be the same as the con- 
figuration of the star with viscosity T2 at time (T2/Ti)i. 
By varying Tvis over a wide range and corroborating this 
scaling, we are confident that the physical behavior we 
observe is real. We can then scale the result of numer- 
ical simulations to the appropriate strength of viscosity, 
provided that the physically relevant viscous timescale is 
much shorter than all the other secular timescales (e.g., 

r ; ' bar / ' 

As discussed in the previous subsection, turbulent vis- 
cosity is likely to dominate the secular evolution. We 
adopt the stress tensor as in Eq. ||7J) and specify the vis- 
cosity rj suitable for modeling turbulence. We consider 
the turbulent viscosity described in [ilj : 

r] ~ p/turbWturb • (36) 

(This viscosity law is also used in some accretion-disk 
models 0.) Typically, Wturb is proportional to the sound 
speed Cs- Hence 77 ~ kurhP\/P/p ~ {knrh/cs)?. For 
simplicity, we assume that /turb / Cs is constant inside the 
star. Then we have 

Tl = ypP . (37) 

where i^p is a constant, and is related to the coefficient of 
kinematic viscosity vhy vp = {p^/ P)v. In our numerical 
simulations we specify the value of vp for each run. 

We model the initial stars as rotating polytropes with 
polytropic index n = 1, so that P = kpq. It is convenient 
to rescale all quantities with respect to k. Since k^/^ ^^s 
dimensions of length, we can define the following nondi- 
mensional variables [43 





= K-^I-'X^ , 


n = 




(38) 


M 




R = 




(39) 


Po 




Dp = 




(40) 
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where the spacetime coordinates are = {t,x,y,z). 
However, to simpHfy our notation, we will drop all the 
bars. Hereafter, all variables are understood to be in 
"k = 1 units." 

Using Eq. (|^ . we can see that Tyis scales with R, p, 
and vp as 



V VpP 



which for n = 1 becomes 

R^ 
pvp 

For definiteness, we take TvIs to be 



Tvis — A- 



PO.maxi^P 



(41) 



(42) 



(43) 



where i?eq is the equatorial radius, po.max is the maxi- 
mum value of pq in the star, and A is a dimensionless 
constant. We use the constant A to approximately match 
Tvis to the rate of decay of differential rotation observed 
by our simulations. With the appropriate Tvis, the value 
of (7^1/ cr^", which is proportional to the rate of energy 
dissipation [see Eq. (|16() ]. is expected to decay like 



(fT^.^T^")|t=oexp(-2Vrvi.) . (44) 



We measure TvIs by numerically evolving a given star, and 
observing the decay with time of < a'^^Up^v >, the average 
value of (J^^CTfj^u throughout the star, weighted by rest 
density. We determine the value of A by requiring that 
Tvis roughly corresponds to the e-folding time of the decay 
of \/<a^^''o7l^7^. Carrying out this measurement on a 
sampling of the stars used below, we find that A « 0.23 
in all cases. Thus we use A = 0.23 to define Tvis in the 
sections below. 

We note that the turbulent viscosity adopted above 
for our numerical treatment is roughly equivalent to 
an "a" model provided we identify a ~ (cslR)vp ~ 
{M/R^Y/'^vp. Eq. for Tvis is then equivalent to 

Eq. jSH). 



F. Radiative cooling 

Our stars do not evolve isentropically. Viscosity heats 
the matter on a timescale Tvis, as shown by Eq. 
At the same time, neutrino radiation carries away heat, 
cooling the star on a timescale TcooI. We will carry out 
simulations below in two opposite limits, which we de- 
scribe in some detail in Appendix ^ In the no-cooling 
limit, Tcooi S> Tvis, so we ignore radiative cooling and 
simply evolve Eqs. (|15|I - H17() . In the rapid- cooling limit, 
T"cooi ^ Tvis, so we evolve Eqs. lfT5)l - lfT7j) as before, but 
without including the viscous heating term in the energy 
equation [Eq. (|16|l ] . This will allow net heating by adia- 
batic compression but not by viscosity. The viscous heat 



is assumed to be (instantaneously) lost by radiation in 
this limit, while viscous braking proceeds. The emitted 
radiation will carry off some momentum as well as en- 
ergy, causing a modification of Eq. IjlTfl . but this will 
have a much smaller effect provided the luminosity does 
not exceed the (neutrino) Eddington luminosity. Baum- 
garte and Shapiro "i^l investigated the loss of angular 
momentum in binary neutron star merger remnants due 
to radiation, and they found it to be fairly small. We 
therefore feel justified in ignoring radiative corrections 
to Eq. l(T7|l. 



G. Numerical Implementation 

1. 2-1-1 Dimensional Code 

Our hydrodynamical scheme employs Van Leer-type 
advection with artificial viscosity to handle shocks. We 
also use a "no-atmosphere" approach, in which the den- 
sity at any point on our grid can fall exactly to zero. 
Our hydrodynamical algorithms are described in detail 
in Paper I. We have evolved the above equations both 
in two dimensions, assuming axisymmetry, and in three 
dimensions. Using axisymmetry saves us computational 
time and allows us to use higher resolution. However, 3D 
runs must still be carried out for rapidly rotating systems 
in order to check for the occurrence of nonaxisymmetric 
instabilities. There are several ways to evolve in axisym- 
metry. One could write the field and hydrodynamic evo- 
lution equations in cylindrical coordinates (zu,z,(p) and 
evolve in this coordinate system. This has the advan- 
tage that one can explicitly remove the dependence of 
the variables on (p. Unfortunately, there are singularities 
in the cylindrical coordinate system which can make the 
evolution of the field equations in these coordinates dif- 
ficult. Instead, we choose to evolve the metric variables 
(7y, Ay, 0, r', a, and /?*) in axisymmetry using 
the Cartoon method [23| . In this approach, variables are 
evolved on a Cartesian grid consisting of three planes cor- 
responding to y — —Ay, y = 0, and y = AY. Then the 
middle (y = 0) plane is evolved using the 3D evolution 
equations in Cartesian coordinates. Each time the middle 
plane is evolved forward one timestep, the y = — AY and 
y = AY planes are updated by applying the assumption 
of axisymmetry. Thus, the value of a tensor f at location 
(tu, z, ±ip) on y — ±AY is equal to f at (w, z, 0) rotated 
by (coordinate) tensor transformation about the z-axis 
by angles ±ip. Since an arbitrary point (tu, z, 0) will gen- 
erally not coincide with any gridpoint on the y = plane, 
interpolation in x is necessary to apply this update. We 
use third-order polynomial interpolation, so that we do 
not lose second-order accuracy. 

With the hydrodynamic evolution equations, we also 
have the choice of either evolving in cylindrical coordi- 
nates or evolving in Cartesian coordinates using the Car- 
toon prescription. Like Shibata in his work on axisym- 
metric star collapse [2^ , we choose to evolve the fluid 
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variables in cylindrical coordinates, i.e. we use Eqs. H22() - 
(|25|l . This is superior to using the Cartoon method be- 
cause Eq. (|22|l can be finite differenced in such a way 
that the total rest mass will be exactly conserved (except 
for flow beyond the outer boundaries). In the absence 
of viscosity, angular momentum also becomes a numer- 
ically conserved quantity. We have found that evolving 
the fluid variables in cylindrical coordinates gives signif- 
icantly more accurate runs than evolving via Cartoon 
hydrodynamics. The drawback of using 2D evolutions 
is the instability caused by the coordinate singularity on 
the a; = axis. (This instability is also present if we 
use the 3D Cartesian Navier-Stokes equations together 
with the Cartoon boundary condition.) There are sev- 
eral ways of removing this instability. Shibata |2^ adds 
a small artificial shear viscosity 



dtCSA/p.) 



(45) 



where A is the flat-space Laplacian, and Si/ p^, = huj 
is the momentum variable evolved by the code of |28| 
(instead of Si). 

We have confirmed that adding such a term to Eq. H17|) 
can stabilize our code. However, since we will be studying 
the effects of real shear viscosity, we instead choose to 
remove the instability using a higher-order dissipation 
scheme. Namely, we add a small Kreiss-Oliger dissipation 
term lil 



dtSA^ 



i6Ar 



(46) 



We use Cko — 0.2 for all the simulations reported in this 
paper. 

In both 2D and 3D simulations, we assume that our 
system preserves equatorial symmetry across the z — 
plane, and we therefore only evolve the z > portion of 
the grid. In 3D runs, we make the added assumption of 
TT-symmetry, which allows us to evolve only half of the 
remaining grid, which we choose to be the y > half. 
When performing 2D runs, we evolve only the region x > 
since the values of variables in the a; < region can be 
deduced from the values on the a; > region from the 
assumed axisymmetry. 



2. Finite differencing 

We compute all spacial derivatives using standard cen- 
tered differencing. We integrate forward in time using a 
3-step iterated Crank-Nicholson scheme. So, for exam- 
ple, we evolve the equation dtf ~ f{f) from timestep n 
at time t to timestep n + 1 at time t+AT by the following 
algorithm: 

(i) Predict : = /" + ^Tf{P) 

(ii) 1st Correct: 

2p+i = pi _^ AT[0.4/(/") + 0.6/(1/"+^)] 

(iii) 2nd Correct: 

jn+l ^ fn ^ Ar[0.4/(/") + 0.6/(2/"+')] 



As discussed in Paper I, the coefficients 0.4 and 0.6 were 
chosen to improve stability. 

In the presence of viscosity, the time-differencing is 
not entirely straightforward, due to the presence of time 
derivatives of in ct^^ , and of time derivatives of cr^j/ in 
the Navier-Stokes equations. Thus, dtSk (which gives 
dfUk) is an expression which itself contains dtUk and 
dfuk- Since the viscosity is a small perturbing force on 
the fluid motion, we find that it is sufficient to split off the 
viscous terms and integrate them separately (operator 
splitting). In particular, we compute dtUk and dta^ ap- 
pearing in the viscous terms in a non-time centered way. 
Consider the S evolution equation for the viscous piece: 

dtS = S{a, a), where we have suppressed all indices. To 
evolve this equation from timestep n to timestep n + 1, 
we need to know the time derivatives of u and a. When 
performing the predictor step, these time derivatives are 
approximated by subtracting values of the fields on the 
timestep n from those on the previous timestep, n — 1. 
(i) Before predictor step, 

"-i]/Ar, 



compute u" — [u" ~ u" 



1/2. 



Note that these time derivatives are centered at n 
We then carry out the predictor step. 

(ii) Predict : ^5"+^ = ^" + AT^(cr", cr"-i/2) 
From the predicted values of u and a, we construct time 
derivatives centered at n -I- 1/2 and use these in the cor- 
rector step. 

(iii) compute from 



-1/2 _ 



l^n+l/2 



u"]/Ar, 

1 =a(iu"+i,iu"+i/2), 



(iv) 1st Correct: 

2 



2^n+l/2 



AT[0.4S'(cr",(T"-i/2) 
-)-0.6^(V"+i 

2 gn+1 

1 -u"]/Ar 

1 = a(2u"+i 2 



(iii) compute ^^"+1 from 

2^«+l/2 ^ r2y 



u- ' 1/2), 
^ (T-- -cr"]/Af 
(iv) 2nd Correct: 

^n+l ^ gn ^ Ar[0.4S'(cr", cr"-l/2) 

+0.6^(2fT"+l,2a"+l/2)] 

By differencing the equations in this way, the dom- 
inant nonviscous terms in the evolution equations are 
accurate to second order (except for small effects due 
to the use of the coefficients 0.4 and 0.6 in the corrector 
steps) , but small viscous terms involving time derivatives 
are only accurate to first order in AT. Numerical con- 
vergence tests show that our code is nearly second-order 
convergent in space and time. We find this to be suffi- 
cient for our purposes. When computing a^^, we first use 
Eq. ||SJ) to get the spatial components . The remaining 
components do/i a-re then obtained from the conditions 
M^cr^i/ = 0. 

As in most other Eulerian hydrodynamic codes, high 
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velocities can easily develop in the low-density regions 
near the surfaces of our stars. The method for evolving 
such regions in the absence of viscosity is described in 
Paper I. Since calculation of the shear tensor involves 
taking derivatives of the velocity field, we are unable to 
calculate it accurately in the very low-density regions. 
To ensure stability, we set = in regions where po < 
10- Po,inax- Since these low-density regions contain an 
insignificant amount of rest mass, this prescription should 
not affect our evolutions. We confirm this by varying the 
cutoff density in several test problems and checking that 
the effect is negligible. 



H. Diagnostics 

Our most important diagnostics are the total mass- 
energy M and angular momentum J. These are both 
defined by surface integrals at infinity |45|: 

M = V77™7^"(7™«,,-7jn,™)d2S. (47) 



7 - — ' 



(48) 



Using Gauss' Law, these surface integrals can be con- 
verted to volume integrals: 



M 



1 



1 ~ 
16^ 



247r 

1 - e'* ~ 



K^) (49) 



167r 



(50) 



Utt ' 167r ' ' / 



In axis ym metry, the angular momentum integral simpli- 
fies to 



J. 



xSyd X 



(51) 



Henceforth, we drop the subscript z since all angular mo- 
mentum is in the z-direction. 

The mass M and angular momentum J in our grid 
should be strictly conserved only in the absence of radia- 
tion (although gravitational radiation carries no angular 
momentum in axisymmetry). Since the energy and an- 
gular momentum emitted in gravitational waves are neg- 
ligible in our runs, this means that M and J should be 
conserved in the no-cooling limit. They thus serve as use- 
ful code checks in this limit. In the rapid-cooling limit, 
the mass computed by the volume integral over the 
numerical grid will not be conserved — thermal energy is 
carried off of the grid by thermal radiation. The expected 
rate of mass-energy decrease due to thermal energy loss 
can be computed by differentiating Eq. (|49|) with respect 



to time. To lowest order, we can ignore the effects of the 
quasi-stationary loss of thermal energy on the spacetime, 
and so ignore the time derivatives of field variables. Then 
only the first term in Eq. I|49|) will be important. 



dM 
~dt' 



cooling 



dt ' ^ 



dt 



(52) 



coolinp; 



where cooling is the component of the time deriva- 

tive of p caused by loss of internal energy due to cooling. 
This quantity may be computed by applying the chain 
rule to Eq. l(T^ : 



dM 




f dp 




de 






~dr 


cooling 


Jv de 








dt 



z^'^d^x 



cooling 



(53) 

Changes in u'^ are ignored because they represent a 
higher-order influence on dM/dt. The rate of change 
in e* due to cooling is given by the effective balance of 
heating and cooling that characterizes the rapid-cooling 
limit (see Appendix A). Thus, 9te^|cooiing is minus the 
value of dtCi, due to viscous heating, i.e. 



de^ 
dt 



cooling 



ae'Mpoe)^'-^'>/^a'^^<j^0 



(54) 



[see Eq. (|16|) ]. From equations (|53|l and H54() . it is 

straightforward to construct 



dM 

~dr 



cooling 



Q/3 



(55) 



r-t- 1 



Finally, the quantity which is nearly conserved in the 
rapid-cooling limit (up to losses due to gravitational ra- 
diation) is 



M + M, 



cooling 



M 



, / dM 
dt' — 
dt' 



(56) 



cooling 



In both the no-cooling and rapid-cooling runs, we can 
divide M into its constituent pieces: the rest mass Mq, 
internal energy mass Mi, kinetic energy T, and gravita- 
tional potential energy W, defined by 147|| 



(57) 
(58) 



(59) 
(60) 

where dV = au^e^'^d^x is the proper 3-volume ele- 
ment. To study the effects of heating, it is useful to 



Mo = 


Jv 




M, = 


1 (poe)dV 

JV 




T = 


il7T°(u°)- 


'dV 


W 


M - Mo - M, 
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break up the internal energy e into its "cold" compo- 
nent eo — Pq~^/{T — 1), and its "thermal" component 
Cheat = e ^ eo- Then we can break up Mi into cold and 
hot components 



M,e = / (poeo)dV 
Jv 

Mih = / (poehoat)dV . 
Jv 



(61) 
(62) 



Note that in the rapid-cooling limit, = 0. 

Finally, we also compute the circulation along closed 
curves. For a closed curve c with tangent vector A'^, the 
circulation is defined to be 



(63) 



C(c) = 6 hu^X^'dC , 



where C parameterizes points on c [i.e. = (d/dC)^]- 
According to the Kelvin-Helmholtz theorem, the circu- 
lation C will be conserved in the absence of viscosity 
if c moves with the fluid and if the fluid is barotropic 
[P = P{pq)]. When viscosity is present or the equation 
of state is more general, as in the case of nonisentropic 
flow, the Navier-Stokes equations give 



d7 



-=-(iA^p„ 



P, 



dc 



(64) 



If 1] = 0, the second term in the integrand vanishes, and 
if P = P{pq), the remaining term is an exact differential. 
Then dC/dr integrates to zero, in accord with the Kelvin- 
Helmoltz theorem. 

In axisymmetry, we choose to evaluate C on circular 
rings on the equatorial z = plane, so that ( = (p. The 
ring c intersects our 2D grid at a point on the a;-axis. By 
our symmetries, the curve c will always remain circular 
and always remain at z — 0, so the Lagrangian point 
representing c only moves in x. Evolution in t and in r 
are simply related by d/dr — u^d/dt. Since the system is 
axisymmetric, the integrand, being a scalar, is constant 
along c, so Eqs. (|53|l and (|64|) simplify to 

C = 2iihu^ = 2TThxUy (axisym) (65) 
dC _ 1 dC _ Att 
It " " T^^'^"'' '''' 

= [(x^ae^Vy"),* + (^'ae«V/),A] (66) 



Hence the quantity Ctot given by 



Ctot — C + Cvis — C — I dt 



dt' 



(67) 



is conserved, even in the presence of viscosity. 

Finally, we compute the Hamiltonian and momentum 
constraint violations [given by equations (16) and (17) of 
Paper I] . We monitor the L2 norm of the violation of each 



constraint. We compute the L2 norm of a gridfunction 
by summing over every gridpoint i: 



(68) 



The constraint violations are normalized as described in 
Paper I [Eqs. (59) and (60)]. 



III. CODE TESTS 

In Paper I, we presented our relativistic hydrodynamic 
code. This code evolves the coupled Einstein field rel- 
ativistic hydrodynamic system on 3D grids, assuming 
perfect-fluid hydrodynamics. We demonstrated the abil- 
ity of our code to distinguish stable from unstable rel- 
ativistic polytropes, to accurately follow gravitational 
collapse of rotating stars, and to accurately evolve bi- 
nary polytropes in quasi-stationary circular orbits. When 
black holes appear on our grid, we can employ excision to 
remove the spacetime singularities from our grid. Tests 
of our black hole excision algorithm using this code were 
reported in (48(| for single rotating black holes in vacuum 
spacetimes and in for black holes that arise during 
the collapse of hydrodynamic matter. In this section, we 
test the adaptations of this code which force axisymmet- 
ric evolution and the modifications which allow a physical 
viscosity. Simulations performed with our axisymmetric 
code show that stable and unstable TOV stars are cor- 
rectly distinguished. The code also achieves approximate 
second order convergence in the evolution of linear grav- 
itational waves and TOV stars. Below, we describe test 
runs on rotating stars in some detail. First, we consider 
stable and unstable uniformly rotating stars, as well as 
a stable differentially rotating star, in axisymmetry and 
without viscosity. We then test the sensitivity of our 
code to nonaxisymmetric dynamical bar formation. Fi- 
nally, we check that physical viscosity is implemented 
correctly by considering stable uniformly and differen- 
tially rotating models. A summary of the models used 
for these code tests is given in Table I. The models are 
initially n = 1 equilibrium polytropes and are evolved 
using a F-law equation of state with F = 2 [see Eq. @]. 
Initial data for all of these models are obtained from the 
relativistic, rotating star equilibrium code of 47j. Stars 
A and B were also studied in Paper I [49j . For the differ- 
entially rotating stars C, D, and E, we choose the initial 
rotation profile 



U Uu 



(69) 



where is the angular velocity of the fluid, Oc is the value 
of f2 at the center and all along the rotation axis, i?oq 
is the equatorial coordinate radius, and the parameter 
A, which measures the degree of differential rotation, is 
chosen to be unity. In the Newtonian limit, Eq. H69|) 
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TABLE I: Initial equilibrium models for code tests. 



Star 


M 


Mo 




a 

PO,max 


T/\W\ 






7? 


A 


0.170 


0.186 


4.10 


0.241 


0.032 


1.0 


155 


0.88 


B 


0.171 


0.187 


3.48 


0.363 


0.031 


1.0 


125 


0.87 


C 


0.183 


0.200 


4.53 


0.155 


0.095 


0.346 


60.6 


0.73 


D 


0.241 


0.260 


5.47 


0.061 


0.234 


0.383 


52.4 


0.37 


E 


0.259 


0.277 


5.92 


0.045 


0.263 


0.381 


57.4 


0.28 



^ Maximum rest-mass density. This does not correspond to the center of the 

star for hypermassive, toroidal models D and E. 

^ Ratio of Q at the equatorial surface to Q at the center. 

Initial central rotation period. 

Ratio of polar to equatorial coordinate radius. 



reduces to the so-called "j-constant" law |5C 



fir 



1 



(70) 



We note that the maximum mass, n ~ I TOV polytrope 
has mass M = 0.164 and compaction Rcq/M = 3.59 [47| . 
All of the axisymmetric tests in this section were per- 
formed with a modest resolution of 64 x 64 and outer 
boundaries at about 12M. Passing these tests success- 
fully with modest resolution helps establish the robust- 
ness of our code. 



A. Tests in axisymmetry 

To demonstrate that our "axisymmetrized" code can 
distinguish stable and unstable uniformly rotating mod- 
els, we consider stars A and B. These stars lie along a se- 
quence of constant angular momentum, uniformly rotat- 
ing stars. As described in Paper I, star A lies to the left 
of the turning point on the M-pc equilibrium curve, while 
star B lies to the rigiit. Then by the theorem of Friedman, 
Ipser, and Sorkin |5l|, star B is secularly unstable to ra- 
dial perturbations, while star A is stable. Since the onset 
of dynamical radial instability is very close to the onset 
of secular instability for such sequences j20| |. we expect 
that star A will be stable to collapse, while star B will 
be dynamically unstable. When evolved in axisymmetry, 
star A persists for more than 7Prot without significant 
changes in structure, where the central rotation period is 



27r 



^c{t - 0) 



(71) 



The oscillations in pc, which correspond to radial pul- 
sations, have an amplitude of < 7%. For this run, the 
Hamiltonian and momentum constraints are satisfied to 
within 2%, while M is conserved to better than 1%. 
Meanwhile, the unstable uniformly rotating star (star 
B) collapses, with an apparent horizon first appearing 
at time t ~ 2Prot, corresponding to 21.3 light crossing 
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FIG. 1: Axisymmetric evolution of uniformly rotating stars. 
Star A (solid lines) is stable, while star B (dotted lines) is 
unstable to collapse. The upper window shows the central 
density normalized to its initial value, while the lower gives 
the central lapse. The solid dot indicates the first appearance 
of an apparent horizon during the collapse of star B. 



times of the grid. At this time, the constraints are satis- 
fied to within 6% and M is conserved to within 3%. Thus, 
stable and unstable uniformly rotating stars are clearly 
distinguished even at this moderate resolution. Figure Q 
summarizes the results for these two runs. 



Next we consider the evolution of a differentially rotat- 
ing star using our axisymmetric code. We evolve star C 
for a time ^ 15Prot- Throughout the simulation, all con- 
straints are satisfied to better than 4.5%, while M, J, and 
Afo are conserved to within 3.5% (J and Mq decrease due 
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to flow beyond the outer boundaries). In the absence of a 
dissipative mechanism to brake the differential rotation, 
the structure of the equilibrium star should not change. 
We find that we can numerically hold this equilibrium 
state for ISPiot- Note that the small amount of Krciss- 
Oliger dissipation employed for numerical stability does 
not alter the rotation profile of the star. After this time, 
inaccuracies at the center, manifested by growth in 
and high-frequency oscillations in po, begin to grow. We 
monitor the evolution of the circulation for three differ- 
ent fluid elements chosen at the following initial locations 
in the equatorial plane: (a) r = i?eq/4, (b) r = i?cq/2, 
and (c) r = 3i?eq/4, where i?oq is the initial radius of the 
star. We find that, for t < 22.5Prot, the circulation is 
conserved to better than 5% for all three of these points. 
After this time, the same inaccuracies cause the circula- 
tion to deviate from its initial value. 

Because viscosity tends to smooth irregularities in the 
velocity field, the problems near the axis and the inac- 
curacy of the azimuthal velocity can be controlled by a 
small shear viscosity. To test this we evolve star C for 
~ 60Prot with a very small shear viscosity (such that 
Tvis — 550Prot)- Because TvIs is so much greater than 
the length of the simulation, the small viscosity does not 
significantly alter the structure of the star. We find that 
the behavior of Q improves considerably, while the small- 
scale variations in po near the axis do not occur. In ad- 
dition, the circulation values for the same three points 
which we studied in the previous case are conserved to 
within 5% for more than 50Prot- Thus, even a tiny shear 
viscosity significantly lengthens the period during which 
our runs are accurate. As we will describe below, the 
presence of a small shear viscosity allows us to evolve ax- 
isymmetric models accurately for hundreds of Prot , cor- 
responding to thousands of M. 



B. Tests of dynamical bar mode sensitivity 

We now demonstrate that our 3D code is sensitive 
to the nonaxisymmetric dynamical bar-mode instability. 
This sensitivity is important because the results of ax- 
isymmetric runs for a particular case will only be valid 
physically if it can be demonstrated that nonaxisymmet- 
ric modes do not develop in the corresponding 3D evolu- 
tion. We consider two models, D and E, which we expect 
to be dynamically stable and unstable to bars, respec- 
tively. This expectation is based on earlier 3+1 fully 
relativistic evolutions of these stars by Shibata, Baum- 
garte, and Shapiro |53 |. who studied the formation of 
bars. Both models are hypermassive, toroidal configura- 
tions with high values of T/\W\ (0.230 for D and 0.258 
for E). These models are identical to models Dl and D2 
considered in j52j |. To test for bars, we add a nonaxisym- 
metric density perturbation of the following form to the 



axisymmetric initial data: 

P = Po + , (72) 

where Sb parameterizes the strength of the initial bar 
deformation. We choose 6b = 0.1 for both models D 
and E. We then re-solve the constraint equations as in 
[53I I to ensure that they are satisfied on the initial time 
slice. The growth of a bar is indicated by the quadrupole 
diagnostic (see jl^l), 

Mo J + 

We will take \Q\ = \JQ*Q as a measure of the magnitude 
of the bar deformation. 

We evolve star D for a time 8.9Prot, during which M 
and J were conserved to within 0.7% and all constraints 
were satisfied to within 2%. For Star E, the run was ter- 
minated after 6.3Prot and M and J were conserved to 
within 1.0%, while constraint violations were < 5.5%. 
Both runs were performed in 7r-symmetry on uniform 
grids with resolution 128 x 64 x 32. The outer boundaries 
in the x-y plane were at 16. GA/ for star D and 19.3Af for 
star E. The results are shown in Fig. [21 This test clearly 
shows the growth of the bar mode for star E, while star 
D does not form a bar even with the substantial initial 
perturbation. 



C. Tests with viscosity 

As a first test of our shear viscosity implementation, we 
demonstrate that uniformly rotating configurations are 
unaffected by the presence of even a large viscosity. We 
evolve the uniformly rotating, stable star A with = 0.2 
(such that Tvis ^ O.OQProt) for -lOOProt - 15, 500M. The 
mass M is conserved to within 0.1%, and J to within 
1.5%. (Note that, even in axisymmetry, J is not iden- 
tically conserved by our finite differencing scheme when 
viscosity is present.) All constraints are satisfied to bet- 
ter than 1.1% for the duration of this run. The result- 
ing evolution of the central rest-mass density and central 
angular velocity are shown in Fig. |3| These quantities 
oscillate on the radial oscillation timescale (~tff) with 
amplitudes of several percent, and this run encompassed 
roughly 120 oscillation periods. Because the oscillations 
are radial, they are not appreciably damped by the shear 
viscosity. For the entire run, the average values of pc and 
drop by about 4.5% and 6.5%, respectively. These 
small deviations are due to the accumulated numerical 
error of the finite differencing and are reduced by increas- 
ing resolution. Since the star does not change apprecia- 
bly over many viscous timescales, our code simulates the 
correct physical behavior for this case. 
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FIG. 2: Quadrupole diagnostic for evolutions of rapidly rotat- 
ing hypermassive stars. Star D is stable and star E is unstable 
to dynamical bar formation. The initial m — 2 perturbation 
decays for star D (dotted line) but grows for star E (solid 
line). Note that, for each curve, the time axis is normalized 
by Prot which differs for the two stars. 



FIG. 3: Evolution of the stable, uniformly rotating star A 
with high viscosity. The central density (shown in the upper 
window) and central angular velocity (lower window) oscillate 
without changing appreciably for over 100 rotation periods 
(15,500A/). 



We now test the viscous evolution of the differentially 
rotating star C. We choose viscosity lyp = 0.015, so that 
Eq. (|43|l gives Tvis ~ 5.5Prot- In axisymmetry, we ran 
this case for 84.5Prot = 15.4rvis = 5, 120M, during which 
time M and J are conserved to within 0.4% while all of 
the constraints are satisfied to better than 1.1%. Figure^ 
shows several snapshots of the angular velocity profile 
in the equatorial plane for the 2D case taken at various 
times. This clearly shows that the presence of viscosity 
drives the star toward uniform rotation. As a quantita- 
tive test of the action of viscosity, we check that the circu- 
lation evolves according to Eq. Ht)6|) . Choosing three fluid 
elements in the initial configuration, we track these fluid 
elements and calculate the circulation C for each one, as 
well as the time-integrated contributions from viscosity 
Cvis [see Eq. (|H7|) ]. The fluid elements are chosen at the 
same locations as for the inviscid test of star C in Section 
ImH The resuhs are shown in Fig. [Sj which gives the 
circulation, the viscous contribution, and their sum, Ctot- 
For all three cases, Ctot is conserved to better than 2% for 
the entire run. Thus, angular momentum is transported 
correctly for many tens of rotation periods (thousands of 
M) when a significant shear viscosity is present. 
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FIG. 4: Angular velocity profiles in the equatorial frame at 
selected times for star C with TvIs ~ 5.5Prot. The presence 
of viscosity drives the star to uniform rotation on a viscous 
timescale. 



We also used this case to test the scaling behavior of 
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FIG. 5: Evolution of the circulation for three selected fluid el- 
ements of star C in axisymmetry (using vp = 0.015). The dot- 
ted line gives the circulation C, the dashed line gives the time 
integrated contribution due to viscosity Cvm, and the solid line 
gives their sum C, which is well-conserved [see Eq.(IS3]. Each 
quantity is normalized by the corresponding initial circula- 
tion. 



FIG. 6: Energy dissipation rate (cr^^cr^y), normalized to its 
initial value, for several runs with star C. The upper win- 
dow demonstrates the proper scaling of our solutions with 
vp. Note that the time axis is scaled according to the appro- 
priate value of vp. The lower window compares {(J^"a^v) for 
runs in axisymmetry (2D) and 7r-symmetry (3D), both with 
i/p = 0.015 = Vfj. These evolutions agree fairly well. 



our solutions with i^p. Results are shown in the upper 
window of Fig. which gives the evolution of {(J^^ff^^^) 
for several values of z/p versus scaled time. We define the 
energy dissipation rate via shear viscosity, {a^'^a^v), as 
in Section III El [see Eq. This quantity decays due 

to the action of viscosity. The figure shows that our so- 
lutions obey the proper scaling with v-p, i.e. they evolve 
identically but on a timescale inversely proportional to 
the adopted viscosity {yp). Hence our results can all be 
scaled to the much smaller viscosities likely to be appro- 
priate for physically realistic viscosity in stars. We pro- 
vide further demonstration of scaling in Section IIV B 21 
(see Fig.ini. 



The results of the axisymmetric run with — vq agree 
fairly well with a 3D, 7r-symmetric run performed for the 
same model. This 3D run employed 64 x 32 x 32 grid 
zones, giving only half of the resolution of the axisym- 
metric run. We terminated the 3D run after 6.35Prot- 
M is conserved to within 0.4%, while the constraint vi- 
olations are < 3.6%. Because of the lower resolution in 
this case, 3.2% of the total angular momentum is lost 
(as opposed to 0.4% for the much longer axisymmetric 
run). A comparison of (cr'^'^iT^i/) for the 2D and 3D cases 
is plotted in the lower window of Fig. |^ and shows good 
agreement. 



IV. DYNAMICAL EVOLUTIONS 

A. Introduction and discussion of models 

Having shown simulations for several test models, we 
now present the evolution of five differentially rotat- 
ing, dynamically stable stellar models in which viscos- 
ity changes the structure of the stars in nontrivial ways. 
Our models are summarized in Table |nl and Fig. We 
first perform short, 3D simulations without viscosity on 
all the five models to make sure that they are all dy- 
namically stable to bar formation. Each of them is then 
evolved with our axisymmetric code in both the rapid 
and no-cooling limits described in Section III Fl The ini- 
tial data for the five stars are again computed with the 
relativistic equilibrium code of [43 ■ The stars obey an 
n = 1 polytropic equation of state P = p^. We adopt the 
rotation law given by Eq. (|69|l with A — 1. This rota- 
tion law has been found to be a good approximation to 
the angular velocity profile of proto-neutron stars formed 
from core collapse 55] . In the case of a binary neutron 
star merger, the remnant can form a dynamically stable 
hypermassive neutron star provided the remnant mass 
does not exceed about 1.7 times the maximum mass of 
a nonrotating spherical star Our adopted rotation 

law is also found to be a reasonably good approxima- 
tion to the angular velocity profile of these hypermassive 
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TABLE II: Initial Models 



Case Mo/Mo,TOv'' 


Mo/Mo,sup" 


M 


R,JM J/M'' 


T/\W\ P,ot/M 




I 1.69 


1.38 


0.279 


4.48 


1.0 


0.249 


0.33 


38.4 


0.2 


II 1.39 


1.13 


0.228 


4.40 


0.85 


0.188 


0.32 


41.3 


0.07 


III 1.39 


1.13 


0.232 


5.54 


1.0 


0.224 


0.37 


54.2 


0.15 


IV 1.39 


1.13 


0.234 


6.27 


1.1 


0.244 


0.31 


63.3 


0.2 


V 1.0 


0.81 


0.168 


8.12 


1.0 


0.181 


0.40 


103 


0.15 



^ If this ratio is greater than unity, the star's mass exceeds the TOV limit for n — 1 polytropes 
(Mo.TOV = 0.180). 

^ If this ratio is greater than unity, the star's mass exceeds the uniformly rotating (supramassive) upper limit 
(Mo, sup ~ 0.221) and is therefore hypermassive. 

" The values of i/p are chosen such that the viscous timescale Tvia ~ 3Prot ~ lOrpp. 



neutron stars 0. 

In all of our axisymmetric calculations, we use a grid 
size 128 X 128 with an outer boundary at 14M for the 
most massive and compact star (star I), and 24M for the 
least massive and compact star (star V). Initially, the 
equatorial radii of the stars are only about Rcq ~ 5M. 
However, viscosity causes the outer layers to expand and, 
in some cases, we find that a few percent of rest mass 
is lost due to material flowing out of the grid. In each 
model, we choose the value of the viscosity coefficient 
Vp such that the viscous timescale defined by Eq. H43|) 
is Tvis ~ 3Prot ~ lOrpF- With this moderate strength of 
viscosity, we need to evolve the stars for 100-200Prot in 
most cases to follow the complete secular evolution and 
determine the final fate of the stars. The reason is that 
in most cases, viscosity generates a low-density envelope 
around the central core. Since our viscosity law has rj cx 
P, the viscosity in the low-density region is small. (The 
density throughout the envelope is greater, however, than 
the cutoff density below which rj = 0.) Hence the effective 
viscous timescale increases with time and it takes longer 
for the stars to achieve a final state. 

Four of the five stars we consider are hypermassive, and 
we expect viscosity to change their structure significantly. 
Star I is the most hypermassive star (Mq = 1.38Afo.sup, 
where Mo^sup — 0.221 is the mass limit for uniformly 
rotating n — 1 polytropes, i.e. for stars at the mass- 
shedding limit 113). We find that this star eventually 
collapses to a black hole, but a substantial amount of rest 
mass is leftover to form a massive accretion disk. We 
consider three other hypermassive models (stars II, HI 
and IV) to study whether or not all hypermassive neutron 
stars will collapse in the presence of viscosity. Stars II, 
HI and IV have the same rest mass {Mq = l.lSA/p.sup) 
which is slightly smaller than that of star I but differ- 
ent angular momenta J. We find that stars III and IV 
never collapse, but evolve in a quasi-stationary manner 
to a uniformly rotating core surrounded by a low-density, 
disk-like envelope. Star II eventually collapses to a black 
hole if we impose rapid cooling. In the no-cooling limit, 
however, this model forms a uniformly rotating core sur- 
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FIG. 7: Rest mass Mo and spin parameter J/M^ for the five 
selected models in Table UTI The dashed line denotes the mass 
limit of uniformly rotating supramassive n = 1 polytropes. 
Mo — Mo, sup. All stars above this line are hypermassive and 
require differential rotation to be in hydrostatic equilibrium. 



rounded by a substantial disk. Star V is the only non- 
hypermassive model. As expected, this star does not 
collapse under the action of viscosity. However, viscos- 
ity cannot drive the whole star to rigid rotation, because 
the angular momentum of the star exceeds the maximum 
angular momentum allowable for a rigidly-rotating star 
having the same rest mass. Instead, viscosity again leads 
to a uniformly rotating core and a differentially rotating 
disk-like envelope. The final outcomes of the five models 
are summarized in Table iHll 

We also performed 3D simulations on stars I and IV to 
search for unstable, nonaxisymmetric secular modes. A 
nonaxisymmetric bar instability usually develops when a 
star is rotating rapidly, i.e. has a sufficiently large T/|VF|. 
Of the five models we study, stars I and IV have the high- 
est r/|IF|. We do not find any nonaxisymmetric insta- 
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TABLE III: Summary of simulations. 



Case 2D/3D Cooling 


ifinal / -Prot ^ 


Initial T/\W\ 


Final T/\W\ 




Fate 




Afo,disk/Mo 


i^disk/ J 


I 


2D 


no 








o. y 


RH -1- disk 


fi 
u.u 


n 9"^ 

u.zo 


fi^ 
u.uo 




2D 


yes 


1 Q 7 




U. 10 


1 n 

l.U 


ijri -i- QISK 


U.o 


n 91 


u.oo 




3D 


yes 


Q 

y .o 






l.U 




u.o 


n 1 8 




II 


2D 


no 


oof? 


A 1 n 
0.19 


U.09 


2.0 


star + disk 




A 1 

O.lo 


U.OD 




2D 


yes 


57.7 




O.W 


1.0 


BH + disk 


0.7 


0.10 


0.36 


III 


2D 


no 


105 


0.22 


0.09 


4.3 


star + disk 




0.21 


0.68 




2D 


yes 


315 




0.12 


1.0 


star + disk 




0.15 


0.58 


IV 


2D 


no 


99 


0.23 


0.10 


7.3 


star + disk 




0.25 


0.76 




2D 


yes 


235 




0.13 


1.0 


star + disk 




0.17 


0.62 




3D 


yes 


11.5 






1.0 


no bar 








V 


2D 


no 


105 


0.18 


0.09 


3.7 


star + disk 




0.13 


0.52 




2D 


yes 


171 




0.13 


1.0 


star + disk 




0.09 


0.38 



The time at which the simulation was terminated. 
^ This quantity corresponds to an average of P/ over the final configuration of the star weighted by rest-mass 
density at the end of the simulation. Thermal pressure generated by viscous heating causes P/ > 1 (recall 
K = 1). We find that the viscous heating is much more significant in the low-density region than in the core. 

These values are obtained by solving Eqs. 17411 - 1821 . 

The quantity r/|W| is undefined when the star undergoes a dynamical collapse. The number given here is an 
approximate value before the star becomes dynamically unstable. 

bilities in these two models, and the 3D results roughly 
agree with the axisymmetric results. 

We discuss the results of our simulations in detail in 
the following subsections. 



B. Star I 



1. 2D with no cooling 

Star I is the most massive neutron star we study. We 
first perform an axisymmetric calculation with no cool- 
ing. At t — 0, the star has a toroidal density profile, 
i.e., the maximum density occurs off center (see the up- 
per left panel of Fig. O . As viscosity gradually brakes 
differential rotation, the star readjusts to a monotonic 
density profile. Figure |H1 shows the maximum rest-mass 
density and the minimum value of the lapse as a function 
of time. Figure ini shows the meridional rest-mass density 
contours at various times. We see that a meridional cur- 
rent is built up in the process. However, the magnitude 
of the meridional velocity (< 0.01c) is much smaller than 
the typical rotational velocity (~ 0.3c). 

Viscosity destroys differential rotation and transfers 
angular momentum to the outer layers. In the early 
phase of the evolution, the core contracts and the outer 
layers expand in a quasi-stationary manner. As the core 
becomes more and more rigidly-rotating, it approaches 
instability because the star is hypermassive and can- 
not support a massive rigidly-rotating core. At time 
t « 27Prot ~ ll^Vis, the star becomes dynamically un- 
stable and collapses. An apparent horizon appears at 
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FIG. 8: Maximum value of rest-mass density (upper panel) 
and minimum value of lapse (lower panel) as a function of 
time for star I in the presence of viscosity. The solid (dashed) 
curves represent the case without (with) cooling. In both 
cases, the central core collapses to a black hole, and leaves 
behind a massive accretion disk. 



time t w 28.8Prot- Without black hole excision, the code 
crashes about lOM after the horizon appears because of 
grid stretching. About 30% of rest mass remains out- 
side the apparent horizon at this point. We then con- 
tinue the evolution using the excision technique described 
in 25:]. We are able to extend the evolution reliably for 
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FIG. 9: Meridional rest-mass density contours and velocity field at various times for star I. The simulation was performed by 
assuming that the system is axisymmetric and experiences no cooling. The levels of the contours (from inward to outward) are 
po/po,max = iQ-° '^^W+°-^) ^ where j = 0, 1, ... 12. In the lower right panel {t = 28.8Prot), the thick curve denotes the apparent 
horizon. 
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FIG. 10: Evolution of the black hole mass Mh and the rest 
mass of the disk Mo, disk after the appearance of the apparent 
horizon at i = 1106M = 28.8Prot. Note that time is plotted in 
units of M (IProt = 38. 4M). Black hole excision is employed 
to track this late evolution. 



another 55M. The system settles down to a black hole 
surrounded by a massive ambient disk. The rest mass 
Mo^disk and angular momentum Jdisk of the disk can be 
calculated by integrating the rest-mass and angular mo- 
mentum density over the volume outside the apparent 
horizon [c.f. Eqs. H57|l and (|51|l ]. The mass of the black 
hole can be estimated by the proper circumference of the 
horizon in the equatorial plane: = Cu/'^t^ (assuming 
that the spacetime can be described by a Kerr metric). 
Figure ^1 shows the evolution of Mh and the rest mass 
of the disk Mo^disk- We estimate the final values of M/,,, 
Afo,disk, and Jdisk by fitting these curves to analytic func- 
tions of the form A+B exp(— Ct) and extrapolating these 
fitting functions to i ^ oo [Sg- We estimate that the 
mass of the final black hole is Mh « 0.82M. The asymp- 
totic rest mass and angular momentum of the ambient 
disk are found to be Mo, disk ~ 0.23 Mq and Jdisk ~ 0.65 J. 
We can infer from the conservation of angular momen- 
tum that the final angular momentum of the black hole is 
Jh « 0.35 J. Hence we find A/M^ ~ 0.52( J/M^) w 0.52. 

The formation of a massive disk is mainly due to the 
fact that viscosity transports angular momentum from 
the inner core to the outer layers. The material in the 
outer region is unable to fall into the black hole because 
of the centrifugal barrier. The final mass of the black 
hole and disk can also be estimated independently from 
the conservation of specific angular momentum using a 
method developed by Shapiro and Shibata |57j , which we 
apply below. 

During the dynamical collapse, the effect of viscosity 
is negligible. Since the spacetime is axisymmetric, the 



specific angular momentum j = hu^ of a fiuid particle 
is conserved. For a Kerr black hole of mass Mh and 
angular momentum Jh, the specific angular momentum 
of a particle at the innermost stable circular orbit (ISCO) 
jisco is given by 



Jisco 



VMhrnUrL - 2aVMhr^s + Q^) 
r^»(r^, - 3Af;,r„s + 2a VM^)!/^ 



(74) 



' msV ms 

where a = Jh/Mh- The ISCO radius is 



r„,s = Mh [3 + Z2- v/(3-^i)(3 + Zi + 2Z2)] , (75) 
where 



Zi = 1- 
and 



1 - 



A. 



1/3 



1 



Jh_ 

Ml 



1/3 



1 - 



Ml 



1/2 



1/3- 



(76) 



(77) 



The rest mass and angular momentum of the escaping 
matter in the envelope with j > jisco is given by 



AJo,disk 
J disk 



]>nsco 



p^i d X 



(78) 
(79) 



J>JISCO 



We assume that the energy radiated by gravitational 
waves is negligible so that the total mass-energy of the 
system is approximately conserved. Hence we have 



M = Mft + Mdisk 

J = Jh + Jdisk ■ 



(80) 
(81) 



For a bound system, the contribution to the mass-energy 
of the matter in the disk, Mdisk, is smaller than its rest 
mass Mq disk- We write 



M, 



disk 



qMo. 



disk 



(82) 



We consider two opposite limits for q: q — 1 and 
q — M/Mq. The value 5 « 1 is a good approximation in 
the weak gravity regime. In the limit where Mh <C M, we 
have Mdisk ~ M and Afp.disk ~ Mq. Hence in this limit, 
q w M/Mo, which is 0.92 for star I (see Table HJ). We 
expect that the correct q lies somewhere between these 
two extremes, which are not very different. The mass 
and angular momentum of the black hole can be esti- 
mated by solving the system of transcendental Eqs. H74|l - 
(|82|l at a particular time slice during the pre-excision dy- 
namical collapse phase, including the time slice at the 
onset of dynamical collapse (where q is close to unity). 
We find that the values of Mh and Jh are insensitive 
to q. They are also insensitive to which time slice we 
choose to do the calculation. We find Mh ~ 0.75M and 
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2. 2D with rapid cooling 
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FIG. 11: Angular velocity profiles in the equatorial plane at 
various times during the evolution of star I. 



Jh ~ 0.35 J {Jfi/M'^ « 0.6). The rest mass in the am- 
bient disk is found to be Mo,disk ~ 0.23Afo. It should 
be noted that this calculation is based on the assump- 
tions that the spacetime around the disk can be approxi- 
mated by a Kerr metric, and that the disk is moderately 
thin and lies in the equatorial plane of the hole. This 
approximation is not reliable when the disk is massive 
(Afdisk ~ -M) • In our case, we find this calculation agrees 
rather well with the actual asymptotic values determined 
by the dynamical simulation with excision. Henceforth 
we will use Eqs. (|71 )l - l|5^ to estimate Mh, Jh, Mo,disk 
and Jdisk whenever an apparent horizon forms. 

Since viscosity is small in the low-density region, it 
takes longer to remove the differential rotation in the 
outer layers. Figure 1111 shows the angular velocity pro- 
files at various times. We see that by the time the inner 
core collapses, the material in the outer layers is still 
differentially rotating. After the dynamical collapse, vis- 
cosity will cause some of the remaining material to slowly 
accrete onto the black hole. 

We monitor the conserved quantities and the con- 
straints during the entire evolution. Since our finite- 
difference scheme preserves the rest mass, the variation 
of Mq can only come from material flowing out of the 
grid. We flnd that Mq is conserved to 0.01%, and angu- 
lar momentum is conserved to 0.1%. The Hamiltonian 
constraint is violated by less than 0.3% before the dynam- 
ical collapse occurs. It increases to 3% by the time an 
apparent horizon appears. The momentum constraints 
are violated by less than 1% before the dynamical col- 
lapse occurs, and increase to 6% by the time an apparent 
horizon appears. 



We next perform an axisymmetric simulation of star I 
with rapid cooling. The dashed lines in Fig. |S1 show the 
time evolution of the maximum rest-mass density and the 
minimum value of the lapse. As in the no-cooling case, 
the inner core contracts and the outer layers expand in 
a quasi-stationary manner. The core then collapses dy- 
namically to a black hole and leaves behind a massive ac- 
cretion disk. Since there is no viscous heating, the whole 
process occurs more quickly than in the no-cooling case. 
The dynamical collapse occurs at time t w 12Pi.ot ~ StvIs 
and the apparent horizon appears at t « 13.5Prot- Fig- 
ure 1121 shows the meridional rest-mass density contours 
at various times. We estimate, by solving Eqs. iTTIIl - llS^ 
at t « 13.5Prot with q = I, that the mass and angular 
momentum of the final black hole are Mh ~ 0.75M and 
Jh « 0.45J [Jh/Ml « 0.8). About 20% of rest mass 
escapes capture to form an accretion disk. 

In Section Fill CI we demonstrated that when the vis- 
cous timescale is significantly longer than the dynami- 
cal timescale, the secular rates of change of all physical 
quantities scale inversely with viscosity. The secular evo- 
lution of star I with rapid cooling is short enough for us to 
perform another detailed scaling test. Figure ^1 demon- 
strates this scaling behavior by evolving star I with four 
different strengths of viscosity vp = 0.4, 0.2, 0.1, and 
0.05. (The curves in Fig. |S1 correspond to i/p = 0.2.) 
We see that the scaling holds until dynamical collapse at 
time t « td{vp)- When t > t^, the evolution of the sys- 
tem is no longer driven by viscosity. We therefore expect 
that the collapse is independent of the strength of viscos- 
ity as long as the viscous timescale is much longer than 
the dynamical timescale. In the lower panel of Fig. 1131 
we demonstrate that it is possible to shift the time axes 
(t t — td) so that the four viscosity runs yield the same 
result when t — id ^ 0, which indicates that viscosity is 
insignificant during the dynamical collapse. The values 
of td are determined by requiring that the scaling relation 
td{i^2)/td{i^i) ~ v\lv2 holds, and that the four curves be 
aligned when plotted against the shifted time t — td(yp^. 
We found that td(vp)/P,ot ~ 6.1, 12.0, 24.08 and 47.75 
respectively for i^p = 0.4, 0.2, 0.1, and 0.05. The fact 
that we are able to find td{vp) that satisfies these re- 
quirements validates our physical interpretation of the 
two phases of evolution. 

To better visualize the effects of viscosity, we follow 
the motions of ten selected Lagrangian fluid elements. 
Figure E| shows the worldlines of these particles. We 
choose the particles to be in the equatorial plane of the 
star. Equatorial symmetry implies that the particles will 
remain in the equatorial plane at all times. The position 
of a fluid particle X satisfies the equations 



d u-{t-X{t)) 
dt ^' ut{t;X{t)) 



(83) 



We label the particles by the initial fraction of rest mass 
interior to the cylinder of radius X. We see that the 
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FIG. 12: Same as Fig.|5|but for rapid cooling. 
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FIG. 13: Evolution ol the maximum rest-mass density ol 
star I for various strengths of viscosity, assuming rapid cool- 
ing. Upper panel: the curves coincide when plotted against 
the scaled time prior to dynamical collapse. Lower panel: 
during dynamical collapse, it is possible to shift the time axes 
\t t ~ tdivp)] so that the curves again coincide, which in- 
dicates that viscosity plays an insignificant role during the 
dynamical collapse phase. 




FIG. 14: The worldlines of Lagrangian fluid elements at the 
equator for star 1, assuming rapid cooling. The cylindrical 
coordinate X of the particles at time t = is chosen so that 
the initial fraction of rest mass ra(X) interior to the cylinder 
of radius X is, from left to right, m=0.03, 0.1, 0.2, 0.3, 0.4, 
0.5, 0.6, 0.7, 0.8 and 0.9. The cross in the diagram denotes the 
location of the apparent horizon at the end of the simulation. 
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FIG. 15: Evolution of the mass M of star I in the rapid- 
cooling limit. The mass is not conserved since the thermal 
energy generated by viscous heating is removed. However, 
the sum of the remaining mass and the mass carried away by 
"radiation", McooUng, is approximated conserved. 



particles with the initial mass fraction m < m» ~ 0.8 
move toward the center and ultimately move inside the 
apparent horizon, while those with m <: move away 
from the center and remain outside the apparent horizon. 
This agrees with our estimates of the rest mass of the 
ambient disk. 

Since there is rapid cooling, the mass is not conserved 
because the thermal energy generated by viscous heating 
is removed, as discussed in Section 111 HI However, when 
we account for the mass-energy carried away by thermal 
radiation, AfcooUng, the total mass Mtot = M + AfcooUng 
should be conserved approximately [see Eq. I|56|) ]. Fig- 
ure El shows M and Aftot as a function of time be- 
fore an apparent horizon appears. The total mass is 
well-conserved except near the end of the simulation, 
where the numerical error arising from the grid stretching 
causes a few percent drop in the mass. 

We monitored the violations of the constraints dur- 
ing the evolution. The violation of the Hamiltonian con- 
straint is ~ 0.1% before the dynamical collapse occurs, 
and goes up to 7% at the time when the apparent horizon 
appears. The violations of the momentum constraints are 
also ~0.1% before the dynamical collapse occurs, and 
increase to 8% at the time when the apparent horizon 
appears. 



C. Other models 

The evolution of star 11 due to viscosity is different 
from that of star 1. Although it is still hypermassive, the 
mass of star 11 is smaller than that of star 1. When evolved 
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FIG. 16: Evolution of the central rest-mass density for star II 
with no cooling (solid line) and with rapid cooling (dashed 
line) . 



in the absence of cooling, the star does not collapse to a 
black hole, but forms a rigidly-rotating core with a low- 
density disk-like envelope. When evolved in the rapid- 
cooling limit, however, the star collapses to a black hole. 
Figure [TBI shows the evolution of the central density. 

In the no-cooling case, Star II has not collapsed to 
a black hole by the end of simulation [t = 286Prot = 
STtvis = ll,800Af), but is setthng to a uniformly rotat- 
ing core surrounded by a massive torus. Figure IT7I shows 
the meridional density contours at the beginning and at 
the end of the simulation. Figure ITSl shows the angular 
velocity profiles at various times. Viscosity drives the star 
to a quasi-equilibrium, rigidly-rotating core surrounded 
by a low-density disk. We cannot exclude the possibility 
that some of the outer material will slowly accrete onto 
the uniformly rotating inner core, eventually triggering 
collapse to a black hole. However the star acquires en- 
hanced pressure support against collapse from viscous 
heating (i.e. P/pg > ''^(O) where k(0) = 1). Hence it 
may no longer be hypermassive with respect to this new 
"hot" equation of state, as the simulation suggests. Dur- 
ing the entire simulation, the star loses 1.2% of its rest 
mass and 4.5% of its angular momentum due to material 
flowing out of the grid. Figure [TT?1 shows the L2 norms of 
the Hamiltonian and momentum constraints defined by 
Eqs. (59) and (60) of Paper 1. We see that the violation of 
all the constraints are smaller than 1% during the entire 
evolution of 286P,.ot = 11800 M. 

The values of the ratio of kinetic to gravitational po- 
tential energy, T/|VF|, for all of the stars we studied de- 
crease with time. Figurel^shows the evolution of r/|M^| 
for star II evolved without cooling. Viscosity transforms 
part of the rotational kinetic energy into heat. It also 
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FIG. 17: Meridional rest-mass density contours for star II 
with no cooling. The upper graph shows the contours at t = 
and the lower graph shows the contours at the end of the 
simulation {t — 286Prot = STtvIs). The contours are labeled 
as in Fig. El 



changes the equilibrium configuration of the star signifi- 
cantly, causing a redistribution of various energies. Fig- 
ure 1211 shows the time evolution of various energies de- 
fined in Eqs. H57|l - Ht)2|) . The mass of the system decreases 
by 1.4% due to a small amount of mass flowing out of 
the grid (not visible in the graph). The rotational ki- 
netic energy T decreases slightly. The contraction of the 
core raises the gravitational binding energy \W\, as well 
as the adiabatic part of the internal energy Mic- Vis- 
cous heating generates the thermal energy Mih, which 
prevents the star from undergoing catastrophic collapse. 

In the rapid-cooling case, star II collapses dynamically 
at time t « 57Prot ~ 17.4rvis. An apparent horizon ap- 
pears at i = 57.7Prot- The mass and angular momentum 
of the final black hole are estimated by solving Eqs. H74(l - 
(|H2l: M,, « 0.88M and A « 0.63 J (Jh/Ml « 0.7). 
About 10% of rest mass is left as an accretion disk. 

The situations for stars 111 and IV are similar. The in- 
ner core contracts in a quasi-stationary manner while the 
outer layers expand. Each system evolves into a rigidly- 
rotating core surrounded by a disk-like envelope. The 
stars do not collapse to black holes at the end of the sim- 
ulations whether or not rapid cooling is imposed. Again, 
we do not rule out the possibility that they might col- 
lapse to black holes when enough material accretes onto 
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FIG. 18; Angular velocity profiles at various times in the 
equatorial plane for star II evolved with no cooling. 



FIG. 20: Evolution of T/\W\ for star II evolved without cool- 
ing. 
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FIG. 19: L2 norms of the Hamiltonian constraint and mo- 
mentum constraints for star II evolved without cooling. 



the inner core. 

From Table IlIII we see that at the end of the simu- 
lation, a large amount of angular momentum is trans- 
ported to a massive disk. For stars III and IV, the rest 
mass of the core Mo,corc is smaller than the rest-mass 
limit of a supramassive star Mg.sup- For star II, Mo.corc is 
slightly smaller than Mo.sup in the absence of cooling, but 
is slightly greater than Mg^sup in the rapid-cooling limit. 
For star I, Mo,corc > -^o,sup in both the rapid-cooling 
and no-cooling cases. This suggests that the fate of a 
hypermassive neutron star depends on whether viscosity 
can create a rigidly-rotating core with Mp^coro > Afo,sup, 
in which case it will collapse. Both viscous heating and 
the star's initial angular momentum play an important 



FIG. 21: Evolution of various energies for star II evolved with 
no cooling. 



role in the final outcome. A hypermassive neutron star 
with higher mass and lower angular momentum is prone 
to collapse, whereas viscous heating tends to suppress the 
collapse. 

Finally, we study the effect of viscosity on star V, which 
is non- hypermassive. As expected, the star does not col- 
lapse to a black hole, irrespective of cooling. The star 
eventually evolves into a rigidly-rotating core surrounded 
by a disk. The fact that the star does not simply be- 
come rigidly rotating without shedding mass is due to 
the fact that viscosity conserves Mq and J. For a given 
Mq and equation of state, there is a maximum value of 
angular momentum Jmax(Afo) above which a star can no 
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longer maintain rigid rotation without shedding mass at 
the equator. In the case of star V, it is apparent that 
J > Jmax- Hence viscosity cannot force the whole star 
into rigid rotation. Similar results were found in studies 
of viscous evolution of differentially rotating white dwarfs 
assuming Newtonian gravitation |4(l |. 



D. 3D tests of bar formation 

The results of the axisymmetric runs described above 
will not be physically relevant if the models are secularly 
unstable to bar formation. We evolved stars I and IV in 
3D to check for the formation of bars. These models were 
chosen because, of our five models, they have the highest 
values of T/|V7|. Viscous heating can lead to expansion 
and, hence, a decrease in r/|M^|. Thus, to further in- 
crease the chance of bar formation, we performed these 
runs in the rapid-cooling limit. We superimposed m = 2 
perturbations on the initial data according to Eq. f72|) for 
both stars with magnitude Si, = 0.1, so that \Q\ = 0.014 
initially. Both runs were performed in 7r-symmetry with 
a uniform grid of size 128 x 64 x 32 and outer bound- 
aries in the x-y equatorial plane at 14. 3A/ for star I and 
17. lA/ for star IV. To reduce computational costs, the 
extent of the grid in the z-direction is only half as large 
as in the x- and y-directions. This setup is feasible be- 
cause these models are initially highly flattened due to 
rapid rotation and because their evolution results in an 
expansion which is largely horizontal. For star I, we find 
that IQI decreases in magnitude until the code terminates 
due to the collapse, when |Q| = 0.0015. Before the col- 
lapse, all constraints are satisfied to within 3.5% while 
M and J are conserved to within 3%. For star IV, \Q\ 
also decreases, reaching 0.0023 after 11.5Prot = 3.3Tvis, 
when the simulation is terminated. In this case, the con- 
straints are satisfied to within 6.0% while M and J are 
conserved to within 1.4%. We also find that the rest den- 
sity contours remain nearly axisymmetric throughout the 
evolutions of both stars. These results indicate that both 
stars are stable against secular bar formation on the vis- 
cous timescale. 



V. DISCUSSION AND CONCLUSIONS 

We have simulated the evolution of rapidly rotating 
stars in full general relativity including, for the first time, 
shear viscosity. Our findings indicate that the braking of 
differential rotation in hypermassive stars always leads 
to significant structural changes, and often to delayed 
gravitational collapse. The rest mass, angular momen- 
tum, and thermal energy all play a role in determining 
the final state. We performed axisymmetric numerical 
simulations of five models to study the influence of these 
parameters. In the presence of shear viscosity, the most 
hypermassive model which we studied (star I), collapses 
to a black hole whether we evolve by ignoring cooling. 



or by assuming rapid cooling of the thermal energy gen- 
erated by viscosity. However, the viscous transport of 
angular momentum to the outer layers of the star results 
in mass outflow and the formation of an appreciable disk. 
Next, we considered three hypermassive models (stars II, 
HI, and IV) with the same rest mass Mq, but different 
values of the spin parameter J/M'^. These models have 
smaller Mq than star I, and are therefore less prone to 
collapse. Star II, which has J/M'^ = 0.85, collapses when 
evolved in the rapid-cooling limit, leaving behind a disk. 
But without cooling, this model evolves to a stable, uni- 
formly rotating core with a differentially rotating massive 
disk. The additional thermal pressure support provided 
by viscous heating prevents collapse in this case. 

In contrast, stars III and IV, which have J/M^ = 1.0 
and 1.1, respectively, do not collapse even in the rapid- 
cooling limit. This is sensible because these models have 
a smaller rest mass than star I, but larger angular mo- 
menta than star II. Though the cores of stars HI and IV 
contract, they are prevented by centrifugal support from 
reaching the necessary compaction to become dynami- 
cally unstable. In both cases, we find low-density disks 
surrounding uniformly rotating cores. However, our sim- 
ulations do not rule out the possibility that slow accre- 
tion of the disk material could eventually drive the uni- 
formly rotating cores to collapse. Disk formation also 
occurs for star V, which is differentially rotating but 
non- hypermassive. Since there exist stable, uniformly 
rotating models with the same rest mass, the braking of 
differential rotation in this case does not result in col- 
lapse. However, differentially rotating stars can support 
larger T/|iy| than uniformly rotating stars. In the case 
of star V, there docs not exist a uniformly rotating star 
with the same (high) angular momentum and rest mass, 
so that mass shedding must take place as viscosity drives 
the star to uniform rotation. In the final state, we find a 
rigidly-rotating core surrounded by a low-density, disk. 

Since results obtained from our axisymmetric code are 
physically reliable only for models which are not sub- 
ject to nonaxisymmetric instabilities, we evolved stars I 
and IV in 3D to check for such instabilities. Previous 
studies in Newtonian gravity have found that the secu- 
lar, viscosity-driven bar instability in uniformly rotating 
stars should set in when T/\W\ > 0.14 0,0]. When gen- 
eral relativity is taken into account, the threshold value 
can be somewhat higher Thus, of all of our mod- 
els, stars I and IV have the best chances of developing 
bars since they have the highest T/|iy|. We introduced 
an initial bar-shaped perturbation and ran these cases in 
the rapid-cooling limit. We found that, in both cases, the 
small initial perturbation decays and no bar is formed. 
This is somewhat surprising since r/|VF| is well above 
0.14 in both of these cases. We plan to address this issue 
in a future report. 

For the evolution of each of our five models, we find 
that a massive disk or torus forms in the final state. The 
disk typically carries ^ 20% of the rest mass of the initial 
configuration. Viscosity transports angular momentum 
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from the interior of the star to the more slowly rotat- 
ing exterior. The exterior regions then expand to ac- 
commodate the additional centrifugal force, forming a 
low-density disk. The inner core becomes rigidly rotat- 
ing and, in some cases, undergoes gravitational collapse. 
The disk, however, remains differentially rotating since 
viscosity acts much more slowly in low-density regions. 
For cases in which black holes are formed, the mass of 
the disk may be estimated by integrating the rest-mass 
density for those fluid elements which have specific an- 
gular momentum j greater than the value at the ISCO, 
jisco [see Eq. (I78II ]. The estimates obtained in this way 
agree reasonably well with the results of our numerical 
simulations. Particularly good agreement was found for 
the case of star I with no cooling, for which we were 
able to extend the evolution some 55A/ beyond the first 
appearance of an apparent horizon. The rest mass and 
angular momentum of the disk surrounding the rotating 
black hole could then be calculated directly and agreed 
well with the estimates. We expect that excision tech- 
niques will continue to be crucial in establishing the final 
fate of systems involving matter surrounding black holes. 

In a recent paper, Shibata [IJ numerically simulated 
collapses of marginally stable, supramassive stars. These 
supramassive models were constructed using polytropic 
equations of state with 2/3 < ti < 2 and rotate at the 
mass-shedding hmit with 0.388 < J/M^ < 0.670. Shi- 
bata found that the collapse of these stars results in Kerr 
black holes and that no more than 0.1% of the initial rest 
mass remains outside of the hole. This result is quite 
different from our finding that disks are usually present 
following collapse. However, the initial data for the two 
calculations arc quite different, as well as our inclusion of 
viscosity. The analysis of takes uniformly rotating, 
unstable configurations as initial data and follows their 
dynamical evolution. Our calculations begin with differ- 
entially rotating, stable configurations and follow both 
their secular (viscous) and dynamical evolution. Viscos- 
ity drives our configurations to uniform rotation. We 
find that massive disks usually form as by-products of 
the formation of uniformly rotating cores. This is due 
primarily to the transport of angular momentum from 
the inner to the outer layers. In addition, all of our mod- 
els have 0.85 < J/M^ < 1-1. (Large angular momentum 
is required to generate a hypermassive neutron star in 
equilibrium.) Since this range is higher than that con- 
sidered in 21], our models more naturally produce disks 

All of the phenomena observed in our simulations fol- 
low from the braking of differential rotation in strongly 
relativistic stars. This may be accomplished by viscosity 
as shown here, but magnetic fields are likely to be more 
important. The fate of the hypermassive remnants of bi- 
nary neutron star mergers may crucially depend on these 
effects. The loss of differential rotation support in such a 
remnant may lead to delayed gravitational collapse. This 
collapse could in turn lead to a delayed gravitational wave 
burst following the quasi-periodic inspiral and merger 



signal Our results indicate that if the remnant is 
not sufficiently hypermassive, collapse may not occur, at 
least not until the star cools by radiating away its ther- 
mal energy. Understanding the evolution of such merger 
remnants could aid the interpretation of signals observed 
by ground based gravitational wave detectors, such as 
LIGO, VIRGO, GEO, and TAMA. In addition, short- 
duration GRBs are thought to result from mergers of 
binary neutron stars or neutron star-black hole systems 
p^. Is^ l . In this scenario, the GRB may be powered by 
accretion from a massive torus or disk surrounding a ro- 
tating black hole. We have demonstrated that such disks 
are easily produced during the evolution of hypermassive 
neutron stars. 

The braking of differential rotation may also be impor- 
tant in neutron stars formed in core collapse supernovae. 
Nascent neutron stars are probably characterized by sig- 
nificant differential rotation (see, e.g., "go", "gJ, 162. l6^ 
and references therein) . Conservation of angular momen- 
tum during the collapse is expected to result in a large 
value of r/|M^|. However, uniform rotation can only sup- 
port small values of T/|V7| without shedding mass ( .3^ , 
Chap. 7). Thus, nascent neutron stars from supernovae 
probably rotate differentially. If the induced differential 
rotation is strong enough, hypermassive neutron stars 
can form. Their subsequent evolution and final fate then 
depends on the presence of viscosity or magnetic fields. 
Such considerations may be im por tant for long-duration 
GRBs in the coUapsar model 64]. In this model, the 
GRB is powered by accretion onto the central black hole 
formed through core collapse in a massive star. 

Several interesting astrophysical systems undergo sec- 
ular evolution in strongly gravitating environments. In 
this paper, we have shown that it is possible to study sec- 
ular effects that occur over many dynamical timescales 
using hydrodynamic computations in full general relativ- 
ity. We consider this an important step toward future 
numerical explorations of secular effects in other con- 
texts. In particular, we plan to incorporate MHD into 
our evolution code, as magnetic braking probably acts 
more quickly than viscosity to destroy differential rota- 
tion in many systems, like neutron stars or supermassive 
stars [g^. Our results have also raised the following in- 
teresting question: Under what circumstances are differ- 
entially rotating, compressible neutron stars with high 
T/|Vt^| unstable to nonaxisymmetric modes? We plan to 
address this issue in a future report. 
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APPENDIX A: VISCOUS HEATING AND 
RADIATIVE COOLING 

In this appendix, we describe the thermal properties 
of our configurations. The dissipation of rotational en- 
ergy by viscosity heats the stars, but they may be cooled 
by radiation (e.g. neutrino radiation). The presence of 
radiation contributes a term R""^ to the stress tensor: 
^totai = + -R''"' with T^'•' given by Eq. O. This 
modifies the equations of motion (|13() to 



(Al) 



where is the 4-force density due to radiation (see |6q). 
The specific entropy s of a fluid element with temperature 
T and number density n — po/m changes when there is 
heating and cooling according to 

ds 

nT— = Thcat - A = Tvis + Trad - Arad , (A2) 
CtT 

where r is the proper time along the element's worldline. 
Here, we have separated the contributions from viscosity 
and radiation to the heating rate. The quantity Fvis is the 
viscous heating rate per unit volume, which comes from 
the aapcr"'^ term in Eq. In terms of the thermal 

energy density Utham and viscous timescale Tvis, Fvis is 
roughly 



^therm /"^vis 



(A3) 



In general, a fluid can be heated and cooled by the pres- 
ence of a radiation field. The energy equation becomes 



(A4) 



where the last equality arises from evaluating m^G^ in 
the comoving orthonormal frame of the fluid, where ■u'^ = 
(1,0,0,0). Then 

G° = Trad - Arad = // diydn{K,I, - r^,) , (A5) 



where the integral is evaluated in the comoving frame and 
Kjy, Jy, and rji, are the opacity, intensity, and emissivity 
at frequency v [6^. For applications of interest here, 
radiation mediates net cooling of the viscous-heated fluid. 
Hence, we can set Frad = for simplicity. 



The first law of thermodynamics 
d{e/n) = -Pd{\/n)+Tds , e = pq{1 + e) (A6) 
and Eq. @ give 



nTds = n^d 



F- 1 



F- 1 



-dn 



(A7) 



Here we define the entropy parameter k by P = ku^ , 
where, in general, k — k(s). The form of Aj-ad will depend 
on the details of the neutron star's microphysics, but it 
must have the property that Aj-ad = when k(s) = kq = 
k(s = 0), where kq is the value of k for the unheated fluid 
(i.e. no emission from a zero-entropy fluid) . Accordingly, 
we replace Eq. HA5|I for A^ad by the following illustrative 
form 



A 



rad 



in[K{s) ~ Ko]/rcooi 



(A8) 



where ^ is a constant and TcooI is the radiation timescale. 
Combining equations (|A2|) - (|A8|) . we find 



(A9) 



dn _T -1 ( C/thcrm ^ri[K{s) - Kq] 



dr rJ- 



'^cool 



In the limit TcooI S> TvIs, radiative cooling is unimportant 
and K increases due to viscous heating. We refer to this 
regime as the no-cooling limit. If TcooI Tvis, then the 
first term on the right hand side of Eq. ljA9p may be 
dropped in relation to the second, giving 



dT 



{k - Kq) = - 



e(F-l)(K-Ko) 



'^cool 



(AlO) 



Thus, K is exponentially driven to kq. We refer to this 
regime as the rapid- cooling limit, whereby the thermal en- 
ergy generated by viscosity is radiated immediately and 
does not heat the gas. In effect, Arad = Fvis in this 
limit. In practice we implement this limit by omitting 
the a^j^^a^" in Eq. ((TB)l . Though we consider only these 
two limits in our analysis, we expect that, in reality, heat- 
ing will dominate in some regimes and cooling in others. 
We treat both limiting cases in our simulations. 
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